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Abstract: We present a general method that allows one to decay narrow resonances in Les 
Houches Monte Carlo events in an efficient and accurate way. The procedure preserves both 
spin correlation and finite width effects to a very good accuracy, and is therefore particularly 
suited for the decay of resonances in production events generated at next-to-leading-order 
accuracy. The method is implemented as a generic tool in the MadGraphS framework, 
giving access to a very large set of possible applications. We illustrate the validity of the 
method and the code by applying it to the case of single top and top quark pair production, 
and show its capabilities on the case of top quark pair production in association with a 
Higgs boson. 
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1. Introduction 

The analysis of scattering processes involving heavy resonances at the Large Hadron Col- 
lider (LHC) is of essential importance. One of the reasons is that the Higgs boson couples 
primarily to heavy particles such as the top quark and the weak bosons. Also many theories 
beyond the Standard Model predict the existence of unstable particles at the TeV scale. 
These unstable particles are not observed directly but they initiate a cascade of decays 
(also called a decay branch) ending up with particles that may be observed in the detector. 
Such processes lead to a rich phenomenology because the spin of the heavy resonances 
and the type of coupling to other fields imply non-trivial angular correlations among the 
final-state particles inside a given decay branch (a phenomenon also called decay spin cor- 
relation effects) or among the particles from distinct decay branches (a phenomenon also 
called production spin correlation effects). 

Nowadays virtually all experimental analyses at the LHC rely on tools to simulate 
scattering events and their reconstruction in the detector. As Monte Carlo generators have 
become an increasingly important tool for collider physics, a large effort has been devoted 
to improving their accuracy in particular, by developing parton- level generators based 
on next-to-leading-order (NLO) QCD matrix elements. In this respect, the simulation 
of production events with heavy resonances at next-to-leading-order accuracy in QCD 
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is particularly challenging, because the decay pattern of these resonances may lead to 
signatures characterized by a large particle multiplicity in the final state. Because of the 
intrinsic complexity of next-to-leading-order amplitudes for large multiplicities of external 
particles, the efficiency of the current algorithms imposes some limitations on the range of 
possible applications. In many instances, some approximations must be made, otherwise 
the generation of events cannot be handled in practice. 

In this context, one well-known simplification is the narrow width approximation, 
which delivers a good accuracy in the case of a resonance with a width T much smaller 
than its mass. In this approximation, intermediate resonances are put on their mass shell, 
significantly simplifying the structure of QCD corrections since radiative corrections in 
the production and in the decay do not interfere in the limit F — )• 0. As a result, QCD 
corrections in the production and in the decay can be handled separately. This scheme 
is implemented in Monte Carlo integrators such as MCFM [||, |2|, ^ for specific processes 
involving the top quark. In that implementation, spin correlation effects are retained at 
next-to-leading-order accuracy (in the limit F — )• 0). 

However, as far as Monte Carlo event generators are concerned, the previous approach 
may not be the optimal one. Indeed, in the case of processes with complicated decay 
patterns, the efficiency of the generation of unweighted events (i.e. events with the same 
weight) becomes a serious issue. This problem of efficiency may be overcome by considering 
a stronger assumption, sometimes called the decay chain approximation Q which consists 
of factorizing the squared amplitude into a production factor and a decay factor. Using 
this scheme, the generation of unweighted events can be split into (1) the generation of 
undecayed events in which these resonances appear as on-shell final state particles and 
(2) the decay of the resonances in each production events. Off-shell effects may partly be 
recovered by smearing the invariant mass of each heavy resonance according to a Breit- 
Wigner distribution, and by reshuffling the other momenta in the undecayed event. In the 
simplest implementation of this scheme, each resonance is decayed according to a uniform 
distribution in its rest frame, in which case spin correlation effects are lost. To gain 
accuracy, kinematics of the decay products can be generated randomly according to the 
squared amplitudes for the decay processes In this latter case, decay spin correlation 
effects are partially preserved, while production spin correlation effects are lost as a result 
of the production vs. decay factorization at the squared amplitude level. 

In the search for the optimal approach to generate events with intermediate resonances, 
a compromise must be reached between the efficiency and the accuracy of the approach. 
In this spirit, Frixione, Laenen, Motylinski and Webber Q have proposed a convenient 
procedure (called the FLMW procedure in this paper) for Monte Carlo generators which not 
only preserves the efficiency inherent to the decay chain approximation, but also includes 
nearly all spin correlation effects at next-to-leading-order accuracy. They implemented 
this scheme in MC@NLO |^, ^] and in Powheg Q for specific processes involving the 
production of top quarks and weak bosons. The accuracy of the scheme was demonstrated 
for these processes. 

Since the work of Frixione et al, there has been a tremendous progress in automat- 
ing next-to-leading-order calculations. A new milestone has been reached by promoting 
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Monte Carlo generators into the era of fully automated NLO event generators |jT^, 
This breakthrough has opened various interesting perspectives on hadron collider phe- 
nomenology, by allowing the simulation of a new class of processes at next-to-leading-order 
accuracy. Among these processes, the production of heavy resonances -such as top quark 
pair production in association with a Higgs boson [11, 13]- are directly relevant for the 
ongoing phenomenology at the Large Hadron Collider. Even though these automated NLO 
Monte Carlo generators feature, in principle, no restrictions on complexity of the process 
and particle multiplicity, in practice the CPU cost becomes enormous for high-multiplicity 
final states. Most of the current tools cannot simulate the full production and decay at 
NLO accuracy in a reasonable amount of time; only the generation of undecayed events at 
next-to-leading order is feasible. 

Some frameworks already exist to decay heavy resonances in undecayed events. How- 
ever the existing tools are either limited to specific processes {e.g. the implementation in 
MC@NLO [P), or they do not provide an acceptable accuracy given the current state-of-the- 
art of simulation techniques {e.g. the program BRIDGE Q or the generic decay routines 
in Pythia []1^ |15| or Herwig [|16|). In this work, we provide a generic and accurate al- 
gorithm to decay heavy resonances in undecayed events generated at next-to-leading-order 
accuracy. We revisit the FLMW procedure and we demonstrate that this procedure can 
be fully automated. Special care is taken to handle the off-shell effects and the reshuffling 
of the momenta. We implement this algorithm in the MadGraphS framework |17] and 
we dub this new tool MadSpin. 



2. Spin correlations in the decay 

We start by introducing some notation. The kinematics X of an event Ex can be parametrised 
by a set of independent variables each of them in the range [0, 1]. The func- 

tion function (/> : x — )■ X that maps the variables x onto the kinematics X is usually 
called phase-space mapping. In the corresponding decayed event Ex^y. (with correspond- 
ing kinematics X^k) an additional set of independent variables y = ., . . . is needed to 
parametrise the kinematics of the decay products. In this work, the variables y for each 
resonance decay are always in one-to-one correspondence with the invariant mass of the 
resonance and the angles defining the direction of the decay products in the rest frame of 
the resonance. 

Next-to- leading-order generators typically produce a list of events, where each event 
Ex is characterized by the fully exclusive kinematics X. The FLMW procedure |^] allows 
one to generate the decays of narrow resonances in the event Ex according to distributions 
that retain both decay and production spin correlation effects. For each undecayed event, 
the procedure goes as follows. 

1. The variables y characterizing the decay of the resonances are generated randomly, 
considering -for each decay- a uniform distribution of the decay products in the rest 
frame of the decaying particle, and the decayed event with kinematics X(jk(x, y) built 
upon these variables is reconstructed. 



-3- 



The differential cross section da/dxdy associated with the decayed process is evalu- 
ated with tree- level matrix elements at the decayed event with kinematics Xdk(x,y) 
and an unweighting procedure with respect to the maximum weight is used to decide 
whether this decayed event should be kept or not. More precisely, if Winax(-^) I'ep- 
resents the maximum value of the differential cross section over the range of possible 
decay configurations y's, a random number r is generated uniformly between and 
1, and the decayed event Xdk(x, y) is retained provided that the condition 

'^'^ (Xdk(x,y)) >rt^^axW (2.1) 



dxdy 



is satisfied. Otherwise the procedure to decay the event with kinematics X is restarted 
from step 1 onwards. This ensures that a decayed event Xdk(x, y) is kept with a 
probability that is proportional to the differential cross section evaluated at this 
event. 

It should be emphasized that only tree-level amplitudes are used when evaluating differ- 
ential cross sections in the FLMW procedure. For production events with an extra hard 
radiation, spin correlation effects are included at next-to-leading-order accuracy. For the 
other events, spin correlation effects are included at leading order accuracy, as no infor- 
mation on the one-loop corrections not factorizing the Born is used to generate the decay 
configurations. 

Frixione et al. also recognised that, in the narrow width approximation, the maximum 
weight W^aj^{X) can be expressed as a product of the differential cross section for the un- 
decayed process d(Tno-dk/c^x evaluated at X and a constant factor W^^ that is independent 
of the undecayed event Ex'- 

W^UX) = ^^(X) X W""^ (2.2) 

They provided analytical formulae for the constants VF^'^'s associated with top quark and 
vector boson {W^,Z) decays. 

In the FLMW approach, off-shell effects are recovered by smearing the mass of each 
resonance according to a Breit-Wigner distribution. This requires to reshuffle the momenta 
of the external particles in the undecayed event. The procedure to perform this momentum 
reshuffling is not specified in [^. 

The achievement in this paper is a generalization of the procedure proposed by Frixione 
et al. to arbitrary processes, and the practical implementation in the MadGraphS frame- 
work. This implementation implies that the code can be used for any processes for which 
the matrix elements can be evaluated with MadGraphS, so it applies to any model p!^, p!9|. 

On the conceptual level, generalizing the FLMW procedure requires to solve two prob- 
lems: 

• a procedure to handle off-shell effects and momentum reshuffling in a generic way 
must be established, 
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• the maximum weight that is used for unweighting the events [see Eq. must be 

determined (at least numerically) on an automated basis. 

We present our solutions for these problems in the two next sections, respectively. 
3. Off-shell effects and momentum reshuffling 

When scattering events with the production of heavy narrow resonances are simulated at 
next-to-leading-order accuracy, the width is typically set to zero, as it simplifies the cal- 
culation. However, the off-shellness of these resonances can have an non-negligible impact 
on the kinematics of the events. As an example, in scattering events with the production 
of W bosons, the invariant mass of the decay products of the W boson lies outside the 
window [niw — 10 GeV, inw -\- 10 GeV] with an 6.5% probability (assuming a Breit Wigner 
distribution). Specific analyses may be sensitive to such effects, in which case off-shell 
effects must be incorporated, at least to some degree of accuracy. 

Even though strictly speaking event generation can be factorized into a production 
phase and a decay phase only in the narrow width approximation, prescriptions can be used 
to recover part of the off-shell effects. In the FLMW procedure, the mass of each resonance 
in undecayed events is smeared according to a Breit- Wigner distribution. The generated 
virtualities of the resonances enter explicitly in the set of variables y characterizing the 
decay configuration, so that the unweighting procedure presented in Eq. ( ^ ) is used to 
capture any deviations with respect to a Breit- Wigner distribution when decaying the 
events. We adopt the same strategy in our procedure^. 

The off-shell effects imply that the masses of some external particles in the undecayed 
event with kinematics X are altered, which requires to modify the external momenta in 
such a way that energy and momentum are still conserved. This modification of the mo- 
menta is sometimes called momentum reshuffling. Practical implementations of momentum 
reshuffling rely on ad-hoc prescriptions. However, for a given process, one strategy may 
be better than the others, as it could preserve some of the features of the distributions of 
events as predicted by the hard scattering amplitude associated with the production events. 
This makes also clear the fact that the best approach to reshuffle momenta depends on the 
process under consideration. 

In our algorithm, we pay particular attention to this aspect: momentum reshuffling is 
performed in an optimized way by using diagram-based information of the tree-level scatter- 
ing amplitude associated with the undecayed events. This is achieved by taking advantage 
of the single-diagram-enhanced multichannel integration procedure in MadEvent |pO|| . 
In that procedure, each channel of integration i uses a phase-space mapping, 

</>i : Xi ^ X(xi), (3.1) 

that is linked to a specific Feynman diagram i. The connection between the diagram and 
the channel is that each invariant associated with a propagator appearing in the diagram is 

^It should be mentioned though that we do not consider effects from the partonic density functions when 
calculating the weight associated with a decay configuration. This may lead to a small systematic bias in 
the distribution of 2 — 1 s-channel events with respect to the invariant mass of the s-channel resonance. 
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in one-to-one correspondence with a variable of integration in the channel. The amplitude 
Ai of the diagram associated with integration channel i is used to set the relative weight 
Wi of that channel point-by-point in the phase space: 

|^.|2 

The connection with momentum reshuffling resides in the observation that a phase- 
space mapping i as defined in Eq. (^]^) provides a convenient way to reshuffle the momenta. 
The canonical variables Xi = xj, . . . ,xl characterizing the original event X are first deter- 
mined before smearing the masses of undecayed particles, by applying the inverse mapping 
(f)^^ on the event with kinematics X. Then, after the masses are modified, the kinematics of 
the reshuffled event X is generated by applying the mapping^ (p^ on the canonical numbers 
Xj ((f)* is the same mapping as except for the masses of the undecayed particles, which 
are already determined in the set y). 

For each undecayed event with the production kinematics X, our procedure to decay 
the event goes as follows. 

1. For each channel of integration i associated with the phase space of undecayed events, 
the squared amplitude of corresponding single diagram \Ai\'^ is evaluated at the kine- 
matics X of the production event (before smearing the mass of the undecayed parti- 
cles). 

2. One channel of integration i with mapping {(pi : Xj — )• X) is chosen randomly, with 
the probability to choose channel i equal to Wi given in Eq. ( p.2[ ). The inverse 
function : X — )■ Xj is used to extract the canonical numbers Xj, which will 
remain unchanged when reshuffling the event. 

3. The variables y characterizing the decay are generated randomly, considering -for 
each decay- a uniform distribution of the decay products in the rest frame of the 
decaying particle and a Breit-Wigner distribution for the virtuality of the resonance. 

4. The momenta in the undecayed event are reshuffled by applying the phase-space 
mapping (j)*{y) : Xj — )• X on the canonical numbers Xj extracted at step 2, where (/>*(y) 
is the same mapping as (pi except that the pole mass of each final-state resonance 
has been replaced by the virtual mass generated at step 3 (hence its dependence 
on y). The resulting reshuffled event is denoted X. The momenta Xdk(x, y) of the 
corresponding decayed event are reconstructed from the kinematics X, augmented 
with the subset of the variables y that correspond to the angular distributions of the 
decay products in the rest frames of the resonances. 

5. The differential cross section da/dxdy associated with the decayed process is evalu- 
ated at Xdk(x, y) and an unweighting procedure with respect to the maximum weight 
is used to decide whether this decayed event should be kept or not. With r a random 

^In some rare cases, this is kinematically impossible, we will comment on this problem later on. 



-6- 



number generated according to a flat distribution between and 1, the event is kept 
if 

^ (Xdk(x, y)) > rWr^^X) (3.3) 

where W-caaxiX) is the maximum weight, of which the determination will be discussed 
in the next section. If the condition in ( |3.3D is not satisfied, the procedure to decay 
the event with kinematics X is restarted from step 3 onwards. 

One difficulty in this approach is that the mapping (l)*{y) : Xj — )• X that is used at step 4 
is not always defined. This typically occurs for events in which a pair of unstable particles 
has been generated very close to threshold. If the difference between the invariant mass 
mi2 of these particles and the sum of the pole masses mi + m2 is of the order of the width 
of one of the particles, the sum of the virtual masses may exceed the original invariant 
mass in the event: mi2 < m\ + rrv^. If the invariant mass mi2 is mapped onto one of the 
canonical variables Xj, the mapping (j)*{y) : Xj ^ X is ill-defined. This indicates that the 
whole approach for factorizing the production of undecayed events from the decay cannot 
be used to get the fine details of event distribution very close to threshold. This is no 
surprise, because in this region of phase space the narrow width approximation that was 
used to generate the undecayed events breaks down. In our implementation, if the mapping 
(/)*(y) : Xj — 7- X is not defined, the procedure is restarted from step 3 onwards. This may 
lead to a migration of events very close to a threshold to events just below the threshold. 
This effect is negligible in the limit Ti/rrii — )> 0. 

One main advantage of the procedure outlined above is that if resonance effects are 
present in the production events, those effects are systematically preserved in spite of 
the momentum reshuffling (due to step 2). Indeed, an invariant mass which is resonant 
is systematically mapped onto a variable of integration in all the dominant channels for 
events in the resonance region. Then, by construction, this invariant mass is not affected 
by the reshuffling procedure. 

4. Estimation of the maximum weight 



The unweighting procedure outlined in Section requires an upper bound on the dif- 
ferential cross section, given by the maximum weig lit Wrai,^{X). In the work of g, the 
maximum weight was determined analytically for decays of weak bosons and the top quark 



by calculating the constant W in Eq. (2^) for each of these decay chains. In this work 



we opted for a numerical estimate of the quantity W^^. Such a numerical estimate of the 
maximum weight is better suited to automate the FLMW procedure to any given process, 
in particular for New Physics Models. 

The simplest solution to get a numerical estimate of the maximum weight is to probe 
the phase space of the decay of the first production event with a large number of points N . 
The largest value of the differential cross section for the decayed process normalised by the 
differential cross section for the production process provides a quantity Wf^ (the index 1 
refers to the first production event) that approximates the required upper bound on the 
maximum weight. However, such an estimate fails to account for finite width corrections 
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Figure 1: Calibration of the relevant parameters in the numerical estimate of the maximum weight, 
based on — > events. Left pane: the average and standard deviation of the Wf^ estimate from 
100 production events are shown as a function of the number of phase-space points N . Right pane: 
the parameter ^ extracted to reproduce the true value of the maximum weight is shown as a function 
of the number of production events m (the number of phase-space points N is fixed to 20k). The 
different curves are associated with different sets of production events. 



to Eq. ( |2.2[ ), and it also suffers from large statistical uncertainties (unless the number of 
phase-space points is really large). Our prescription is slightly more elaborate to deal 
with these two issues: 

• we extract the estimates Wf^, . . . , for the first m production events, 

• we express the estimate of an upper bound on W^^, denoted by W^^^, in terms of 
the mean value (Wf^') and the standard deviation std (Wj*^*^) associated with the 
numbers Wf^, . . . , W^: 

<ax = (Wt'') + e std (Wf^) . (4.1) 



In this way, we account for finite width corrections by considering multiple events, since 
those corrections depend on the kinematics of the production event. Statistical uncer- 
tainties as well as fluctuations due to finite width corrections are minimised by using the 
average and standard deviation in Eq. ([1.1|), rather than (for example) the largest of all 
Wf-^. This results in a reliable unweighting efficiency. 

Our prescription contains three adjustable parameters: the number of phase-space 
points N associated with the decay, the number of production events m, and the prefactor 
^ in Eq. ([4.lD . These parameters were calibrated so that the quantity Wmax is as close 
as possible to (while strictly larger than) the true value of the maximum weight, in order 
to optimise the Monte Carlo unweighting efficiency. Based on empirical studies involving 
different processes, we verified that the same calibration can be used for all processes. 

In Figure [l|, we illustrate the calibration of the three parameters A^, m and ^ based on 
pp ^ W events. First, the number of phase-space points is increased until the average 
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of the Wf^ is constant with respect to the number of phase-space points^ (left pane). This 
identifies the number of phase-space points to be used. Next, the true value of W^^ is 
extracted by considering a very large number of production events. Once this value is 
known, the prefactor ^ in Eq. ( |4.1| ) can be extracted for a given number m of production 
events so that the resulting estimate reproduces the true value W^^ (right pane). 

We take the number of production events m large enough so that the associated value of 
^ converges to a constant value. That constant then gives the value of thereby fixing 
the remaining parameter in our prescription. In the code, we have made the conservative 
choice {N =10'^,m = 20, C = 4). 

5. Validation 

In order to support phenomenological studies in an optimal way, Monte Carlo generators 
must show a high level of efficiency, while reflecting as much as possible the accuracy of 
state-of-the-art calculations. Regarding the efficiency, the approach discussed in this paper 
has been implemented in the MadGraph5 framework, and hence it can be used for a 
large class of processes within or beyond the standard model in a complete automated way. 
Moreover, since the decay of a specific event typically requires only a few evaluations of 
tree- level matrix elements, generating the decay is in general fast. The efficiency and the 
flexibility of the tool is at the cost of two approximations inherent to the procedure at 
work: 

1. some finite-width effects may be lost in the distributions of events, as the procedure 
is built upon the narrow width approximation in the first place, 

2. only tree- level matrix-elements are used to calculate the weight associated with a 
specific decay configuration. 

We comment on the validity of these two approximations in this section. 
5.1 Finite widtii effects 

We verified the validity of the first approximation by considering a large class of tree- 
level processes involving heavy resonances. We generated events in two different ways (a) 
we used MadGraphS to generate undecayed events, which were subsequently decayed 
using MadSpin, (b) we used MadGraph5 to generate decayed events, considering a 
finite width for the heavy resonances at all stages of the generation (and including also 
the non-resonant contributions). We compared a large number of distributions involving 
the transverse momentum, the angular separation, and the invariant mass of the final 
state particles. We obtained very good agreement between distributions resulting from the 
generation procedures (a) and (b). Some differences were observed in the tails of resonant 
invariant mass distributions or close to threshold, as expected. 

^While the average becomes constant at large A'^, the Wf'^ do fluctuate event-by-event due to finite width 
corrections. This is displayed by the nonzero standard deviation at large A''. 
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Figure 2: Distribution of events with respect to the invariant mass of the muon pair (left pane) 
and with respect to \/I (right pane) in pp — >■ /i+/i~e+z^e events. The two histograms result from two 
distinct procedures to generate the events: either MadGraph5 is used to generate pp ZW^ (LO 
accuracy) which are subsequently decayed using MadSpin (solid histogram) or MadGraphS is 
used to the generate pp — > fj,'^ fj,~ Ve in one shot, so that all finite-width effects are systematically 
included (dashed histogram). 



To illustrate the largest observed deviations, let us consider the specific example of 
diboson production at the LHC in the channel pp — t- ZW~^ — ?• fi^ fi~e~^i'e. Figure ^ 
shows the distribution of events with respect to the invariant mass of the muon pair (left 
pane) and with respect to -v/I = the invariant mass of the colliding partons (right pane) 
resulting from the previously-mentioned procedures to generate events: (a) MadGraphS is 
used to generate pp — )• ZW^ events (at leading order accuracy) which are subsequently 
decayed using MadSpin (solid histogram), (b) MadGraphS is used to the generate pp — )• 
fi'^ fi~ e'^ Ve events in one shot (including also the non-resonant diagrams with the photon 
splitting 7* — )• p'^ p~), so that all finite-width effects are systematically included (dashed 
histogram). In this last case, we impose the cut m(/i+,/i~) > 40 GeV. 

As expected, procedure (a) fails to reproduce the correct distribution of events far away 
from the resonance region m{fi~^ , fi^) ~ mz, as the distribution of events in these regions is 
sensitive to the non-resonant contributions involving the photon splitting 7* — )• p^p^ ■ We 
observed though that the distribution of events with respect to the invariant mass of the 
muon pair is accurately reproduced in a rather extended region around the pole mass mz- 
although the Z boson is generated on its mass shell in undecayed events, off-shell effects 
are recovered to a very good accuracy when decaying the events. We also observe a good 
agreement for the distribution of events with respect to the invariant mass of the colliding 
partons (\/i), except below the threshold region \/I ~ mz + mw where the effects from the 
finite widths of the Z and W bosons are of primary importance and cannot be reproduced 
in the narrow width approximation. An ad-hoc approach to improve the description below 
the mz + m\Y threshold would be to allow for a (small) change in \/I when the ratio 
{mz + mw)/V§ is close to one. We leave this for future work. 
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Figure 3: Next-to-leading-order cross sections differential in coscj), cos(6'+) and pt{1^) for pp tt 
events in the dileptonic channel. The angles are defined in the text. 

5.2 Spin correlation effects in NLO events 

For the decay of Les Houches Events generated at next-to-leading-order accuracy, the 
procedure in MadSpin retains spin correlation effects at tree- level accuracy, i.e. no in- 
formation from the virtual amplitude is used in calculating the weight of a specific decay 
configuration. The validity of this second approximation can be assessed in some specific 
cases by comparing distributions of events against predictions including spin correlation 
effects at a higher level of accuracy. One possibility is to generate these predictions with 
MCFM (provided that it includes the process at work) which includes all QCD corrections 
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Figure 4: Next-to-leading-order cross sections differential in cos{9) and pt{1^) for single-top pro- 
duction in the four flavor scheme (t-channel only). The angle 9 is defined in the text. 



when calculating the weight of a decay configuration. A limitation of such a comparison 
is that MCFM provides only fixed-order predictions, whereas MadSpin is embedded in a 
scheme involving NLO calculation matched to parton shower. A more advanced validation 
procedure {e.g. involving NLO predictions in the complex mass scheme matched to parton 
shower) is beyond the scope of this paper. 

In this section, we present comparison plots with the results from the Monte Carlo 
program MCFM for tt and single-i production at the LHC at ^/s = 8TeV. For single- 
t we used the four-flavour scheme and considered production via the f-channel. In 
both processes the t and t decay semi-leptonically (b quark -|- lepton -|- neutrino). QCD 
corrections in the decay itself are ignored. The relevant masses are mt = 172.5 GeV and 
rrib = 4.75 GeV. Jets are reconstructed by means of the anti-ZcT algorithm [^], with i?cut = 
0.4 and (pT)min = 25 GeV. For ti we have used the first PDF set from MSTW2008nlo68cl 
(v5.7) [^]. In that case as{Mz) = 0.12018, which evolves by two-loop running. For 
single-i we used MSTW2008nlo68cl_nf4 (v5.8.4) instead ||], where as{Mz) = 0.11490. 
The renormalisation and factorisation scales are set equal to the same scale fi. In the case 
of ti, n is given by the average of the t and i transverse masses, while in the case of single-t 
we use four times the b transverse mass, this choice being based on the work of p5|| . In 
both cases, top quarks are decayed in the semi-leptonic channel. Events are showered and 



hadronised with the program Herwig |16| 



For top quark pair production, the observables that are most sensitive to spin correla- 
tion effects were discussed in 126]: 

• cos{(j)), where (p is the angle between the direction of flight of in the t rest frame 
and the direction of flight of /~ in the t rest frame, 

• cos(0±), where 9+ (9-) is the angle between the direction of flight of (/^) in the t 
(t) rest frame and positive beam direction. 

As can be seen in Figure ^, predictions for these observables using the scheme outlined 
in this paper (aMC@NLO-|-MADSPiN-|-HERWiG) are overall in very good agreement with 
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those generated by MCFM. A good agreement is also found for the pT spectrum of the 
positively-charged lepton [Figure ^, bottom right]. 

For single-top production, the observables that are most sensitive to spin correlation 
effects were discussed in |^^. In the case of t-channel production, the angle 6 is defined 
in the t rest frame as the angle between the directions of flight of and the hardest 
non-&-tagged jet. In this case as well, predictions using the scheme outlined in this pa- 
per (aMC@NLO-|-MADSPiN-|-IlERWiG) are in good agreement with those generated by 
MCFM (see Figure ^). The apparent difference in normalization of the two curves in the 
plot on the left-hand side of Figure |^ is due to the definition of the observable. For the 
aMC@NLO-|-MADSPiN-|-HERWiG predictions there are slightly more events with a non-6- 
tagged jet compared to the fixed order MCFM results (80% versus 76%, respectively). 



6. Application: top-quark pair production in association with a light 
Higgs boson 

In order to illustrate the capabilities of the tool, we apply it to the case of top-quark pair 
production in association with a light Higgs boson at the LHC (running at 8 TeV), con- 
sidering both the scalar and pseudo-scalar hypotheses for the Higgs boson. Due to the 
large irreducible QCD background, any search strategy for this Higgs production process 
relies strongly on the accuracy of the Monte Carlo predictions. QCD correction to these 



processes has been analysed by two groups [11, 13 [ and a comparison between these inde- 



pendent calculations has appeared in Ref. [|2^. In these works it was shown that the NLO 
corrections are very mild, in particular on shapes of distributions. 

To the best of our knowledge, the problem of retaining spin correlation effects in 
events generated at NLO accuracy has not been addressed yet for these processes. This 
problem is trivially solved using the scheme proposed in this paper: NLO parton-level 
events are generated with aMC@NLO, (LHC at 8 TeV, PDF set = MSTW2008(n)lo68cl, 
ruH = mA = 125 GeV, = jip = {rriTiH / A)rnT{t)mT{t))^^^^\ no cuts) and then decayed 
with MadSpin before they are passed to Herwig for shower and hadronisation. In this 
illustration, top and anti-top quarks are decayed semi-leptonically, whereas the Higgs is 
decayed into a pair of b quarks. 

Figure |5| shows the normalized distribution of events with respect to cos((^) (which 
was defined in the previous section), and with respect to the transverse momentum of the 
hardest positively-charged lepton. Although spin correlation effects significantly distort 
the distribution of events with respect to cos (</>), their impact on the pT spectrum of the 
leptons is milder, except at large transverse momentum. The relatively larger effect in the 
tail of this distribution can easily be understood from the fact that the inclusion of the 
spin correlations is a unitary procedure: a small change at low py, where the cross section 
is large, needs to be compensated by a larger (relative) effect at high pT- 

It is interesting that spin correlations have a much more dramatic influence on the 
shape of these distributions than NLO corrections: the leading order results fall directly 
on top of the NLO results for these normalized distributions (both without spin corre- 
lations), as can be seen by comparing the dotted blue and dash-dotted red curves. We 
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COS(l|)) 

Figure 5: Next-to-leading-order cross sections differential in pt{1^) (left pane) and in and cos0 
(right pane) for ttH events with or without spin correlation effects. For comparison, also the leading- 
order results without spin correlation effects are shown. Events were generated with aMC@NLO, 
then decayed with MadSpin, and finally passed to Herwig for shower and hadronisation. 
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Figure 6: Next-to-leading-order cross sections differential in pt{1^) (left pane) and in and cos0 
(right pane) for ttA events with or without spin correlation effects. Events were generated with 
aMC@NLO, then decayed with MadSpin, and finally passed to Herwig for shower and hadroni- 
sation. 



can therefore conclude that preserving spin correlations is more important than including 
NLO corrections for these observables. Naturally, the inclusion of both, as is done here, 
is preferred: it retains the good features of a NLO calculation, i.e. reduced uncertainties 
due to scale dependence (not shown), while keeping the correlations between the top decay 
products. 

The results for the pseudo-scalar Higgs boson are shown in Figure |^ The effects of the 
spin correlations on the transverse momentum of the charged lepton are similar as in the 
case of a scalar Higgs boson: about 10% at small pT, increasing to about 40% at px = 200 
GeV. On the other hand, the cos(0) does not show any significant effect from the spin- 
correlations. Therefore this observable could possibly help in determining the CP nature of 
the Higgs boson, underlining the importance of the inclusion of the spin correlation effects. 
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7. Conclusion 



In this paper, we discussed the decay of events in a LHE file. We showed that the decay 
can be treated in a completely generic way in the MadGraph5 framework. The procedure 
is particularly efficient to generate unweighted decayed events, while preserving spin corre- 
lation effects at tree-level accuracy. Some features associated with the finite width of the 
resonances to be decayed are also successfully reproduced. The practical implementation of 
this procedure -dubbed MadSpin- makes use of the user-friendly MadGraphS interface, 
so that it can be trivially applied to a very large category of processes. 

The most obvious application is that MadSpin provides a very convenient way to 
decay events generated at next-to-leading-order accuracy, as it is naturally embedded in 
a Monte Carlo scheme. Although only tree-level matrix elements are used to unweight 
the decay configurations, we have shown for specific processes that this procedure captures 
essentially all spin correlation effects as predicted by a full next-to-leading-order calculation. 
Therefore such a tool -in combination with a next-to-leading-order Monte Carlo generators- 
is expected to give a practical access to event generation with an improved accuracy in many 
instances. As an illustration we applied the tool to the case of top quark pair production 
in association with a Higgs boson at the LHC. 
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A. Manual 

The MadSpin program is part of the MadGraphS distribution from version 2.0.0 onwards. 
As a consequence, the MadGraphS website can be used for downloading the code but 
also to submit questions or report bugs. The code can either be used in standalone mode 
or called by another program of the MadGraphS Suite (i.e. MadEvent or aMC@NLO). 
In standalone, one launches a MadSpin session by typing: 

. /MadSpin/madspin 

This opens a prompt, similar to that of the MadGraphS Suite. The following commands 
are then available: 

• import FILE: imports the LO or NLO Ihe file and reads the associated header in 
order to load relevant information, e.g. the relevant model. 
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• define LABEL = PARTI PART2 ... : allows to define an additional multi- 
particle tag. The multi-particle tags defined in the header of the file are automatically 
recognized. 

• decay PROCESS: specifies the decay branch initiated by an unstable particle to 
be decayed by MadSpin. Multi-particle tags can be used for the final state particles. 

• set OPTION VALUE: allows to change some internal options of MadSpin such 
as the seed or the value of the maximal weight in the unweighting procedure. Type 
"help options" in the interface for more details on the various available options. 

• launch: runs MadSpin according to the specified options/decay channels. 

• help COMMAND: provides detailed information on a specific command. 

As stated above, MadSpin can also be called directly via the MADEvENT/aMC@NLO script 
(via the launch command or the ./bin/generate_events script). In this case, the script will 
request some information about the programs to run. The exact question depends on the 
programs installed on your computer and on the QCD order of the process (LO or NLO), 
but it typically looks like this: 

Which programs do you want to run? 

/ auto : running existing card 

1 / parton : Madevent 

2 / pythia : MadEvent + Pythia. 

3 / pgs : MadEvent + Pythia + PGS. 

+10 / +madspin: adding MadSpin to the run (before running Pythia) 
[0, auto, 1, parton, 11, parton+madspin, 2, pythia, 12, ... ] 

If you choose to run MadSpin {i.e. if you enter one of the following answers: 11, 12, 13, 
parton+madspin, pythia+madspin, pgs+madspin), the program will propose, in a second 
step, to edit the madspin_card.dat. This file contains the command lines associated with the 
MadSpin program that we just described. Note that MadGraph5 will figure out which 
event file has to be decayed, so that you do not need to include the "import" command 
in this file. After these questions have been answered, MadGraphS proceeds with the 
Monte Carlo generation: (1) it generates the LHE event file before decay, (2) it executes 
MadSpin and (3) if requested, it submits the decayed events to a shower /hadronisation 
program. 
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